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Abstract. We introduce the Statistical Hadronization with Resonances (SHARE) suite of 
programs and perform a study of particle multiplicities as well as bulk thermodynamic quantities 

■ for the RHIC Au-Au reactions at ^/snn = 130 GeV. We also show that the statistical 
f^} hadronization model, with parameters fitted to match pion, proton and hyperon ratios, in 

■ turn correctly and consistently reproduces rapidity particle multiplicity. 

o 

^ e— I 1. Introduction 

The statistical hadronization model Ej has been used extensively to study soft strongly 
interacting particle production since the 1950s. When the full spectrum of strongly interacting 
resonances is included [3], this approach is capable to describe the abundances and spectra of the 
produced particles in detail. The emitted particles carry information about the gross features 
of the hadron source. Their study allows precise extrapolation to unmeasured particles and/or 
kinematic domains, allowing understanding of bulk properties of hadronizing matter such as 
particle multiplicity, mean energy per particle, specific per baryon entropy etc. Once hadron 
spectra are understood, information about the dynamical properties such as collective flow of 
hadronizing matter becomes accessible. 

Other important physical information is contained in the statistical model parameters 
obtained fitting measured particle ratios. These offer additional insights, for example one 
may want to relate the chemical (particle production) freeze-out temperature to the phase 
transformation temperature of the deconfined quark-gluon plasma (QGP) phase into hadrons. 
The chemical freeze-out temperature also helps understand the percentage of total particles 
which originated from resonance decays. The chemical freeze-out parameters fix e.g. the charged 
and neutral hadronic multiplicities, flavor density, baryon stopping, and so on. 

To properly address these issues, standardization of the technical and mathematical tools 
employed in statistical hadronization studies has to occur. Furthermore the model differences 
have to be understood and the different versions have to be unified. Here, we note the chemical 




Figure 1. Chemical non-equilibrium can arise out of kinetic evolution (left), and in rapid QGP 
freeze-out into hadrons (right). 



equilibrium version, where all (light and strange) flavor yields are assumed to have evolved for a 
long time in the hadron phase reaching equilibrium yield |3 El • More refined approaches allow 
for an under-saturation of the strangeness quantum number [7j, implemented quantitatively by 
a parameter 7 S (assumed to be < 1) which affects both s and s in the same way (7^ = 7j, in 
contrast to the fugacity A s = A^ 1 ). The simplest physical scenario applicable here is that the 
hadron source does not live for long enough for strangeness production to approach chemical 
equilibration (Fig. ^ left), while the more rapidly evolving light flavor yields had time to 
equilibrate. 

The third model approach is to allow for chemical non-equilibrium also in the light flavor 
yield. This situation is most likely to be found in a reaction scenario involving formation of the 
QGP. Namely, the phase space density of flavors in QGP is higher than of a equilibrated hadron 
gas in the physical domain explored today at SPS and at RHIC 8 . This is illustrated in Fig. 
^in the right panel, where we show that the phase space density difference compensation leads 
to a step-up in the quark flavor occupancy parameter 7. 

Linked to the question of chemical equilibration is the timescale between hadronization (the 
moment at which degrees of freedom become hadronic) and freeze-out (the moment at which 
hadrons stop interacting). It is generally assumed that there are two freeze-outs, chemical 
and thermal, the former corresponding to particle production and latter to spectral shape 
generated by elastic hadron-hadron scattering. Chemical equilibration requires slow evolution 
of the hadronic gas system, hence to a situation in which hadronization and freeze-out are well 
separated, and considerable modification of all hadronic yields could occur in the interacting 
hadron gas phase {3j. In such a scenario one also expects that the thermal freeze-out is clearly 
different from the chemical freeze-out. 

If QGP hadronizes, the difference between quark and hadron phase spaces makes it likely 
that a sudden freeze-out will occur creating initially hadrons in chemical non-equilibrium; 
In particular, the high entropy density of a QGP phase with massless degrees of freedom 
accompanied by the rapid matter flow driven by the high internal pressure makes it reasonable 
that hadrons are produced in a fast phase transformation, which leads to flavor over-saturated 
state [5], with more qq and ss pairs than is expected at chemical equilibrium. It is necessary to 
ascertain phenomenologically whether, in fact, one needs non-equilibrium to account for particle 



abundances in heavy ion collisions. The Statistical Hadronization with Resonances (SHARE) 
package has been developed allowing to resolve this ambiguity as one of the tasks 

In the following section we give introduce the statistical model, and the details implemented 
in SHARE. We refer the user to Ref. 9 for a detailed presentation and user's manual. This 
is followed in section El by discussion of our results, addressing in turn particle ratios and bulk 
properties. 



2. Main SHARE features 

Many of our experimental friends present at this meeting can obtain on back of an envelope 
the statistical yield of stable hadrons. SHARE extends this to include important refinements 
capable to change these results significantly. 



2.1. Particle yields allow quantum statistics and chemical non- equilibrium 
The Grand-Canonical statistical prescription assumes that enough particles of each flavor are 
produced to keep fluctuations of each quantum number small, and the system volume reduces to 
a "large" normalization constant. This requirement, together with entropy maximization, leads 
to the Fermi-Dirac or Bose-Einstein distribution functions for densities of particle species v. 

n(mi,gi;T,Ti) = ni = gi f : , , (1) 

J (2vr) J T -i exp( x /p 2 + mf/T) ± 1 
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In Eq. (J2J), the upper signs refer to fermions and the lower signs to bosons, respectively. T, is the 
fugacity factor, and mi is the particle mass. The quantity gi = (2Jj + 1) is the spin degeneracy 
factor as we distinguish all particles according to their electrical charge and mass. The index 
i labels different particle species, including hadrons which are stable under strong interactions 
(such as pions, kaons, nucleons or hyperons) and hadron which are unstable (p mesons, A(1232), 
etc.). The second form, Eq. P|). expresses the momentum integrals in terms of the modified 
Bessel function K%. This form is practical in the numerical calculations and is used in the 
SHARE code. Although in principle in some limiting cases this is not a convergent expansion, 
this is a rare exception: the series expansion (sum over n) converges when Tje _mi//T < 1. 
Violation of this condition occurs in practical context only for the pion case within the range of 
parameters of interest. 

In the most general chemical condition 1 , the fugacity is defined through the parameters 
Xji, X q , X s , X c (expressing, respectively, the isospin, light, strange and charm quark fugacity 
factors), and 7 g ,7 s ,7 c (expressing the light, strange and charm quark phase space occupancies, 
= 1 for absolute yield equilibrium). The fugacity Tj is then given by: 

= Xji (x glq ) N « (x s% fi (x clc f° (Vr ff )*S ( W* ( W 3 > (3) 



where 
and 



X q — Xg , Xg — Xg 1 , A c — Ag 1 , (4) 



lq = lq, Is = Is, 7c = IS- (5) 



1 This condition is commonly called chemical non-equilibrium. However, the conventional equilibrium in which 
existent particles are redistributed according to chemical potentials is maintained here. The non-equilibrium 
regarding particle production is, in precise terms, called absolute chemical (non)equilibrium. 



Here, iVj and iV* are the numbers of light (u,d), strange (s) and charm (c) quarks in the 
ith hadron, and NJj, Ng and are the numbers of the corresponding antiquarks in the same 
hadron. 

2.2. Particle yields from resonance decays 

At first, we consider hadronic resonances as if they were particles with a given well defined mass, 
e.g., their decay width is insignificant. All hadronic resonances decay rapidly after freeze-out, 
feeding the stable particle abundances. Moreover, heavy resonances may decay in cascades, 
which are implemented in the algorithm where all decays proceed sequentially from the heaviest 
to lightest particles. As a consequence, the light particles obtain contributions from the heavier 
particles, which have the form 

ni = 6 2 ^i ••• b N -> N -in N , (6) 

where bk^k-i combines the branching ratio for the k — > k — 1 decay (appearing in with 
the appropriate Clebsch-Gordan coefficient. The latter accounts for the isospin symmetry in 
strong decays and allows us to treat separately different charged states of isospin multiplets of 
particles such as nucleons, Deltas, pions, kaons, etc. For example, different isospin multiplet 
member states of A decay according to the following pattern: 

A ++ ^ir++p, (7) 
A + ^(7r + +n) + ^(7r°+p), (8) 

A°^±(n-+p) + ^+n), (9) 
A~^yr~ + n. (10) 

Here, the branching ratio is 1 but the Clebsch-Gordan coefficients introduce another factor 
leading to the effective branching ratios of 1/3 or 2/3, where appropriate. 

To implement this procedure on every particle, one needs to keep in mind that the partial 
widths (product of branching ratio with total width) are often not sufficiently well known. 
In addition, in case of weak decays, an experiment-specific acceptance coefficient is needed to 
correctly model the observed particle rate. This introduces implementation dependent variances 
of statistical hadronization. To combat invisible model variations we suggest an open-source 
"standard" , where resonance decay trees are kept on record and can be updated in a transparent 
fashion, and weak acceptances can be set to correspond to needs of each experiment. This is 
the SHARE code 0. 

In SHARE, as a rule, all decays with the branching ratios smaller than 1% are disregarded. 
In addition, if the decay channels are classified as dominant, large, seen, or possibly seen, the 
most important channel is taken into account. If two or more channels are said to be equally 
important, we take all of them with the same weight. For example /o(980) decays into tttt 
(according to ^U] this is the dominant channel) and KK (according to |1U| this is the seen 
channel). In our approach, according to the rules stated above, we include only the process 
/o(980) — > tttt. Similarly, ao(1450) has three decay channels: wir (seen), 717/ (958) (seen), and 
KK (again seen). In this case, we include all three decay channels with the weight (branching 
ratio) 1/3. A table of allowed decays and branching ratios is provided and can be updated. 
Users can modify this table to study the magnitude of systematic error introduced by incomplete 
knowledge of both resonance masses and decay parameters. 



2.3. Hadron yields allowing for finite resonance width 

If the particle i has a finite width Ti, the thermal yield of the particle is more appropriately 
obtained by weighting Eq. Q over a range of masses to take the mass spread into account: 



hf = / dMn{M,gi\T,Ti) — ^ ^ > n u for Ti -> 0. (11) 

1 7 V ' y ' ' ; 2vr(M-mi) 2 + r2/4 *' 1 V 7 

The use of the Breit-Wigner distribution with energy independent width means that there is a 
finite probability that the resonance would be formed at unrealistically small masses. Since the 
weight involves a thermal distribution n(M,gi] T, Tj) which would contribute in this unphysical 
domain, one must use, in Eq. (|11|). an energy dependent width. 

The dominant energy dependence of the width is due to the decay threshold energy phase 
space factor, dependent on the angular momentum present in the decay. The explicit form can 
be seen in the corresponding reverse production cross sections |11[ I12j. The energy dependent 
partial width in the channel i — > j is to a good approximation: 



,(M) = 6^r 



M 



for M > rriij. (12) 



Here, rriij is the threshold of the decay reaction with branching ratio For example for the 

decay of i := A ++ into j := p + 7r + , we have rriij = m p + m n +, while the branching ratio is 
unity and the angular momentum released in decay is lij = 1. From these partial widths the 
total energy dependent width arises, 

r i -^r i (M)=^r i _ i (M). (13) 

For a resonance with width, we thus have replacing Eq. (jllj) : 
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— V f°° dMn(M, 9i : T,Ti) r^(^) (14) 

Ni^Lv 1 ' (M -mi) 2 + [r i (M)] 2 /4' 



and the factor iV (replacing 2tt) ensures the normalization: 

n-XC^M-^Zmw (15) 

In principle, Eq. (|14|) does not take into account the possibility that the state into which 
one is decaying is itself a unstable state in a thermal bath. Doing this would require a further 
average over the width distribution of receiving state. This higher order effect is at present not 
implemented in SHARE. 



2-4- Canonical effects 

We did not address in SHARE refinement specific to small particle numbers where the grand 
canonical phase space description fails |14| . Thus SHARE is geared to describe particle 
multiplicities where within the causal interaction domain there is an effective 'particle' heat 
bath. To be more specific, the grand canonical phase space will yield correctly e.g. the O yield 
in the limit that the number of strange quark pairs significantly exceeds those required to make 
this particle, e.g. we need about 10 ss pairs. Without doubt these are the conditions prevailing 
in the physical environment we are exploring here, since at RHIC 130GeV about 8 strange quark 
pairs are produced per each participating baryon [El EI!, an d baryon density is about 25 per 
unit of rapidity. 
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Figure 2. Statistical significance profiles for volume for different centralities (left), and it's 
correlation with the charged particle number (right). 

3. RHIC 130 fit results 

3.1. Particle ratios 

We study particle ratios taken at different though similar collision centralities, e,g, 5%, 6%, 
8%, 10%. In first and hopefully good approximation we expect that the physical conditions 
established are similar in all these case. In other words we expect that for particle yields in the 
statistical model to be consistent: 

• The volume should appear as a normalization constant. Temperature and the other thermal 
parameters should be, within error, independent of the centrality for a range of centralities. 

• This normalization constant should present an approximately proportional dependence on 
the total charge multiplicity. 

These requirements can be tested by fitting the total particle multiplicity of a range of 
experimentally measured centrality bins, together with a set of hadron ratios (Fig We 
consider seven most central bins, as measured by STAR, |17| . and find no significant variation 
of the fitted bulk parameters (T, \ q>s , 7 9)S , Xjs) and their errors. The only parameter in the fit 
which does vary from bin to bin is the absolute normalization, needed to describe the charged 
multiplicities; It's statistical significance profiles are shown in Fig. EI (left) ; As noted elsewhere 
(HI j introducing 7 g significantly and consistently improves the statistical significance of the 
fit. 

When the fitted volume normalization parameter is plotted against the mean charged particle 
multiplicity, an approximately linear dependence is observed for all models (Fig. fright). Light 
quark non-equilibrium decreases the necessary volume by as much as 30 %. While this is 
predictable, since the non-equilibrium model invariably leads to quark over-saturation (7^ > 1) 
at high energies and thus higher particle density [S]. The smaller volume required in non- 
equilibrium has not as yet been subject to an in-depth analysis, particularly in light of the 
HBT puzzle |18j : It is sometimes thought that an explosive hadronization scenario is ruled out 
by HBT measurements [Ej- However, as this study the non-equilibrium particle production 
is yielding a higher density than that predicted by equilibrium hydrodynamics models, which 
should rectify this discrepancy. 

We have used RHIC ratios (130 GeV) to calculate thermal parameters within each model. 
While a plot, showing a model-data comparison, is given in Fig. |3 and table Q we refer the 
reader to ^H] for an in-depth data and statistical significance analysis of the fits. 
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Figure 3. The best fits to hadron yields at RHIC-130 GeV. See tabled for parameters of the 
fits and statistical significance. 



Table 1. RHIC-130 GeV hadronization parameters 



parameter 


r = o 


vary 

r pm 


lq-s 

r = o 


= 1 

r ing 


lq = 1, 7. 

r = o 


varies 

r [101 


T [MeV] 
10 4 (A 9 - 1) 
A s 2 


133 ± 10 


135 ± 12 


158 ± 13 


0.157 ± 15 


152 ± 16 


153 ± 23 


708 ± 342 


703 ± 337 


735 ± 390 


730 ± 382 


724 ± 373 


721 ± 363 


1.03132973 


1.03203555 


1.02636207 


1.02788897 


1.0295346 


1.0300848 


7? 


1.66 ± 0.013 


1.65 ± 0.030 


1 


1 


1 


1 


7s 


2.41 ± 0.61 


2.28 ± 0.46 


1 


1 


1.17 ± 0.30 


1.10 ± 0.25 


10 4 (A/3 - 1) 


30 ± 305 


28 ± 293 


59 ± 564 


53 ± 508 


64 ± 525 


59 ± 481 




fit relevance 


N-p = DoF 


16-5 


16-5 


16-3 


16-3 


16-4 


16-4 


X 2 /DoF 


0.4243 


0.4554 


1.0255 


0.8832 


0.6067 


0.7301 


significance 


0.9461 


0.9307 


0.4225 


0.5705 


0.8385 


0.7232 



3.2. Bulk quantities 

We proceed to calculate ^Jsnn = 130 GeV Au-Au bulk quantities. The results are shown in 
table El It is immediately apparent that many of the observed quantities depend strongly on 
whether light quarks are in equilibrium. By contrast, introducing incomplete strange chemical 
equilibration while maintaining jg = 1 does not produce a significant shift in the calculated 
quantities from the full equilibrium values. This is due to the fact that in this case the fit leads 
to practically fully equilibrated system as within error j s = 1 results. 

Overall, the greatest change is that the entropy density is considerably higher within the 
7 9 7^ 1 model. This is consistent with the model's physical motivation, since the high 7 g is 
invoked to conserve entropy in a QGP^ HG transition, without need for a mixed phase and 
expansion. 

Strangeness per entropy, being primarily sensitive to input ratios such as K/tt (fitted), is not 



Table 2. RHIC 130 GeV Au-Au best fit bulk quantities ratios, including charged particles 
(/i ±0 , energy (E), Baryon/antibaryon {B,B), strangeness (s,s) and entropy (S). The bottom 
ratio refers to the fraction of totally emitted particles which arise out of weak decays (assuming 
negligible Kl — ► 37r detectability) 







vary 


7</,s 


= 1 


lq = 1, It 


varies 


ratio 


r = o 


r [To] 


r = o 


r [To] 


r = o 


r ma 


(h+ + h~)/(B-B) 


30.0827 


30.821 


18.613 


20.087 


22.149 


22.758 


(h + - h )/{B - B) 


0.618 


0.618 


0.603 


0.602 


0.632 


0.624 


E/(B -B) 


37.220 


38.152 


31.756 


33.403 


36.004 


36.441 


( s + s)/(B-B) 


9.771 


9.419 


6.328 


6.517 


8.3504 


7.862 


S/(B - B) 


362.679 


369.857 


236.14 


252.001 


280.8769 


284.407 


{h+ + h~ + h°)/(B - B) 


71.265 


70.996 


53.551 


54.479 


62.108 


60.678 


(h + + hr + h°)/(h + +hr) 


1.6427 


1.64 


1.6338 


1.638 


1.6414 


1.641 


{h + -h-)/{h+ + hr) 


0.0142 


0.0143 


0.0184 


0.0181 


0.0167 


0.0169 


E/(h + + h-) 


0.858 


0.882 


0.9689 


1.004 


0.952 


0.986 


{s + s)/{h + + hr) 


0.225 


0.218 


0.193 


0.196 


0.221 


0.213 


o i yit ^ it j 


8 360 


8 546 


7 905 


7 576 


7 493 


7 696 


E/bT 


1.748 


1.796 


1.982 


2.054 


1.943 


2.01 


{B-B)/hr 


0.0469 


0.0471 


0.0624 


0.0612 


0.0534 


0.0552 


(s + s)/h~ 


0.4588 


0.443 


0.395 


0.401 


0.451 


0.435 


S/hr 


17.0297 


17.417 


14.738 


15.493 


15.161 


15.719 


(h+ + h- + h°)/h- 


3.346 


3.341 


3.342 


3.349 


3.352 


3.354 


(B + B)/h- 


0.256 


0.256 


0.327 


0.323 


0.289 


0.295 


S/{s + s) 


37.118 


39.266 


37.314 


38.669 


33.636 


36.176 


S/E 


9.7443 


9.694 


7.436 


7.544 


7.801 


7.804 


(h++h-+h°) weak 
(h++h-+h a ) total 


0.325 


0.326 


0.425 


0.427 


0.410 


0.404 



significantly affected by model choice [161120] . Quantities such as particle multiplicity per baryon 
(which follows entropy per baryon), net charge per baryon, and strangeness per baryon, however, 
vary considerably between models, making them promising probes for statistical model variants 
involving chemical (non)equilibration. We have also show (last line in tabled) that the fraction 
of particles coming from weak decays varies considerably from model to model, decreasing as 
non-equilibrium is introduced. This is to be expected, since all non-equilibrium fits yield over- 
saturated phase space occupancies. These enhance all particles, including resonances, thereby 
increasing the percentage of particles emitted in strong decays. This feature should also provide 
a test for non-equilibrium: Since all weak decays, including those which can not be reconstructed, 
tend to occur at a measurably large distance from the primary vertex, a precise estimate of the 
number of particles which do not come from the primary vertex could also serve as a potential 
probe for equilibrated emission. 

Due to the higher particle multiplicity in the non-equilibrium case, the thermalized energy per 
particle goes down. While this model does not include collective (transverse and longitudinal) 
flow, and hence can not address the total energy in the system (some of which is contained in the 
collective motion, rather than the internal energy), the discrepancy in the energy per particle 
can become a strong constraint if the models under consideration are coupled to a model (such 
as the commonly used "Blast Wave") incorporating flow. Successful fits of abundances and 
spectra within the same model have been made both in the equilibrium and non equilibrium 
cases (2lJ OS [22] , and it will be interesting to rigorously test those for full energy conservation. 



In conclusion, using the SHARE package we have analyzed the Au-Au y / s/v r /v = 130 GeV 
experimental output with a variety of statistical models. All are able to fit particle ratios 
and total charged multiplicity, achieving the highest statistical significance with the full non- 
equilibrium ansatz. We have calculated bulk thermodynamic properties of the system from the 
fitted parameters, and found that many of these are model dependent, with particularly strong 
discrepancies between full non-equilibrium (7g )S 7^ 1) and the rest. We hope that higher statistics 
and 200 GeV data |13j will be able to falsify some of these models unambiguously, and hence 
allow a precise determination of the statistical properties of the system. 
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